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Excitation of surface plasmon waves in extrinsic graphene is studied using a full-wave electro¬ 
magnetic field solver as analysis engine. Particular emphasis is placed on the role played by spatial 
dispersion due to the finite size of the two-dimensional material at the micro-scale. A simple instruc¬ 
tive set up is considered where the near field of a wire antenna is held at sub-micrometric distance 
from a disk-shaped graphene patch. The key-input of the simulation is the graphene conductivity 
tensor at terahertz frequencies, being modeled by the Boltzmann transport equation for the valence 
and conduction electrons at the Dirac points (where a linear wave-vector dependence of the band 
energies is assumed). The conductivity equation is worked out in different levels of approximations, 
based on the relaxation time ansatz with an additional constraint for particle number conservation. 
Both drift and diffusion currents are shown to significantly contribute to the spatially dispersive 
anisotropic features of micro-scale graphene. More generally, spatial dispersion effects are predicted 
to influence not only plasmon propagation free of external sources, but also typical scanning probe 
microscopy configurations. The paper set the focus on plasmon excitation phenomena induced by 
near field probes, being a central issue for the design of optical devices and photonic circuits. 
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I. INTRODUCTION 

One of the main challenges for current device electron¬ 
ics is the full exploitation of the terahertz (THz) electro¬ 
magnetic (EM) spectrum, bridging the gap beitween the 
microwaves and optics. In this context, graphene and 
graphene-derived materials possess a number of unique 
EM properties such as tunable conductivity and slow- 
wave features that may be used to develop high perfor¬ 
mance THz device^®. More generally, low-dimensional 
systems with honeycomb-like geometry and related as¬ 
semblies or hetero-structure are undergoing massive in¬ 
vestigation as materials capable of supporting plasmon 
propagation in the THz frequency rangd^ 12 ! 

A THz beam impinging on a metal atomic force micro¬ 
scope tip has been used to generate guided THz waves 
on graphene in a recent experiment^ the nanometric 
curvature of the tip makes light scatter on the length 
scale of tens of nanometers 2 , providing the EM field to¬ 
gether with the wave-number components required to ex¬ 
cite short wavelength surface plasmons. A similar situa¬ 
tion has been encountered at lower frequencies in scan¬ 
ning microwave microscopy (SMM)P^, where a metallic 
tip can act as a resonant antenna whose tip-termination 
is directly coupled with the sample. 

In this paper, we focus on numerical propagation and 
excitation of plasmon surface waves in spatially disper¬ 
sive graphene as seen from the macroscopic EM perspec¬ 
tive. We use the fact that a momentum matching be¬ 


tween the field source and the surface plasmon can be 
achieved, e.g., by grating or prism couplers. In partic¬ 
ular, we present an EM scheme that ideally recalls the 
set-up of Ref . 13 and take a resonant wire antenna placed 
at sub-micrometric distance from a disk-shaped graphene 
patch. 

Our analysis suggests that the correction to the con¬ 
ductivity response given by spatial dispersion is strongly 
required, at least in excitation phenomena induced by a 
strong spatial gradient of the EM field on the sample. 
We therefore asses that spatial dispersion is an impor¬ 
tant limiting factor to the achievable coupling between 
surface-plasmons and the near field of an SMM tip, be¬ 
cause it produces a broadening of the system-response 
over a wide range of wave-vectors. The high impact of 
spatial dispersion in propagation and modal analysis of 
plasmons has been previously scrutinized^^. However, 
a similar study has not been conducted so far in near 
field excitation problems, which are of central importance 
in numerous practical applications including microscopy 
and photonics. Here, we provide such a characterization. 

From the numerical point of view, we simulate surface 
plasmon excitation in an open computational domain 
avoiding periodic boundary conditions and plane-wave 
excitations 5 . The typical high charge mobility of carbon 
nanostructures leads to a huge increase of the macro¬ 
scopic inductance, i.e., a high reactive energy stored per 
unit length and, accordingly, a slow-wave propagation. 
Consequently, the plasmon velocity can be much smaller 
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than the speed of light with very small wavelengths as 
compared to that of the coupled feeding antenna. This 
also results in a high confinement of the EM field in the 
direction transverse to the propagation, which further 
increases the aspect ratio. We then apply two differ¬ 
ent numerical methods to validate the results, i.e., the 
method of moments (MoM) and the method of finite ele¬ 
ments (FEM). In this way, we provide a comparison of the 
different levels of approximation used in semi-classical 
approaches to the graphene conductivity, based on the 
Boltzmann transport equation for the Dirac electrons and 
the Kubo-formula 17 21 . 


II. THEORY 

The dispersive conductivity of graphene on the THz 
regime is well defined in the literature 115 16 : basically, 
an electric field frequency below a few THz induces a re¬ 
sponse of the valence (7r) and conduction (77*) electrons 
of the material with an energy close to the Fermi energy. 
The dominant part of the response in extrinsic (doped 
or gated) graphene is given by intra-band excitations, 
which can be treated by the Boltzmann transport equa¬ 
tion under the relaxation-time-approximation (RTA) or 
Bhatnagar-Gross-Krook (BGK) model. 

The RTA replaces the relaxation dynamics of each one- 
electron state by a simple exponential decay. Then, the 
transport relaxation time is approximated with the life¬ 
time of the state. Such a formalism provides an accurate 
description of spatial dispersion effects in doped graphene 
with extrinsic Fermi energy shifts below ^ 0.5 eV. Ac¬ 
cordingly, the RTA conductivity is entirely determined by 
the drift currents that arise from one-electron transitions 
and collective modes within the 7r and 7r* bands. 

The BGK model is more general than the RTA, be¬ 
cause it allows for an extra degree of freedom, which 
enforces charge conservation and accounts properly for 
electron diffusiorl 15 ^. 

We begin by specializing to the RTA conductivity 

o-RTA(q) = crt TA (q) + o’raA(q) 

under an applied electric field with momentum q and 
frequency /. cr RTA includes the contributions of elec¬ 
trons (with charge — e and wave-vector k) that occupy 
the 7r*(+) and 7 r(-) bands. The band energies e±(k) are 
populated according to the Fermi Distribution (FD) 

, 1 
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at the absolute temperature T and chemical potential fi. 
We adopt the widely used convention of setting the Dirac 
point energy of graphene, i.e. the intrinsic Fermi level, 


to zero energy. In this way, (i coincides with the Fermi 
level shift caused by the doping (or gating). 

The intra-band components of cr RTA have the 
dyadic (tensor) form 15 


<dk(q) 


Tt r d 2 /: 4w-/±( k ) v ±( k ) 

2tt 2 J Cj — v±(k) • q 
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where: (i) the band energies are linearized as e±(k) = 
±hv F k , with v±(k) = Vk£±(&) representing the elec¬ 
tron velocities and v F = |v±(k)| the Fermi velocity; (ii) 
f f £±(k)-n 1S the energy-derivative of the FD, i.e., 
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(iii) Co = u) — 77 is a complex frequency that includes the 
angular frequency uj = 2 irf and a small shift along the 
imaginary axes, which corresponds to the electron damp¬ 
ing rate 7; the latter depends on the average relaxation 
time r as 7 = 27 r/r. 

In Eq. 0, the (two-dimensional) first Brilluoin 
Zone (BZ) integration is well defined on circular areas of 
the k-space centered at the Dirac points, where the linear 
approximation for the band energies is valid. Neverthe¬ 
less, in most practical uses, an infinite cone-structure is 
assumed for the valence and conduction bands; in other 
words, the wave-vector integral is performed over the 
whole k-space for analytical convenience and a factor 
of 2 is included to account for the inequivalent Dirac 
point. Then, a change of variable from wave-vector 
to the energy leads to the non-dispersive conductivities 

^ta(o) = 4 t ra r in which 
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and I denotes the 2 x 2 identity matrix. The total intra¬ 
band conductivity at q —>• 0 is then 


intra 
r 00 


Sintra — ^intra + 
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( 2 ) 


Interestingly enough, the RTA conductivity 0 and its 
q —>> 0 form Q can be derived from the Kubo formula in 
the optical limit 17 21 . 

Another part of the graphene conductivity is related 
to inter-band processes between the 7r and 7r* bands, and 
is also included in the Kubo formulismPSHUl, though it 
does not contain the dispersive term v± (k) • q. Indeed, it 
has been pointed out that for surface waves supported by 
isolated graphene sheets, and working frequencies below 
a few THz, the spatial dispersion effects on inter-band 
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transitions can be neglecte d 15 * 16 1 . Due to the absence 
of the dispersive term, the inter-band conductivity is a 
scalar, which can be expressed as 


_ ie 2 f°° de / £ _ M - /_ £ _ M 

ainter " nh J 0 to 1-(||) 2 


( 3 ) 


we obtain the tensor components of the intra-band con¬ 
ductivity 


^(q) = Orta ^) 1 + °\ 


= ± 
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which depend on the scalar conductivities 


Therefore, the total conductivity reads 

er(q) = cr RTA (q) + o- inter I. (4) 


In the q —>> 0-limit the latter tends to the non-dispersive 
Kubo conductivit}ff^^ 


&K — ^intra T Winter- (3) 

Now, looking at the denominator in Eq. ([l]), it is clear 
that when v±(k) q and u are comparable, spatial disper¬ 
sion cannot be neglected in surface plasmon excitation. 
This happens particularly in problems where a resonant 
behavior of the EM field is concerned. It is the case of 
some of the examples reported in this work, where the 
interaction between a radiating antenna and a graphene 
patch takes place via near field coupling. 

To give an idea of the numbers involved, the spatial 
harmonics of the excitation field become significant in 
the conductivity response at an operating frequency / of 
the order of ^ 1 THz, for an applied wave-vector value q 
larger than ^ 1 /im -1 . At the same frequency, the elec¬ 
tric field wavelength is about ~ 300 /im. However, the 
near field distribution between the tip of the antenna 
and the graphene sample varies on a sub-micrometric 
scale, depending on the tip radius and distance from the 
sample. Then, the slow-wave effect featured by plasmon 
propagation implies a wavelength reduction of more than 
one order of magnitude, and ensures the matching and 
coupling with the exciting near field. 

Eq. Q can be simplified and made more explicit by 
expressing the planar wave-vectors q and k in polar co¬ 
ordinates, say, q = ( q,0 q ) and k = (&,$/-). Then, as 
shown in the appendix, we can reduce it to the following 
expression 

. ie 2 v 2 r°° 

^(q) = J o dkk fe ± (k)-v ( 6 ) 

cos 2dk sin 2d^ 
sin 2idk — cos 2dk 
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where the BZ integral has been turned to a wave-vector 
integral over the whole k-space (with the factor of 2 from 
the inequivalent Dirac points being included). 

After some straightforward manipulations on the angu¬ 
lar integral in Eq. 0 that are reported in the appendix, 
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Therefore, summing over the d= channels, we find 


0-RTa(<?) = Ct/ta(9) + O-RTA id) = 


^intra 


yi - v 2 g 2 /w 2 


and 


cr RTA (d) = o-+ TA (g) + cr RTA ( q ) 

_ (V 1 -^ 2 /^ 2 - 1 ) 2 gointra 

y/l — v 2 q 2 /uj 2 v F q 2 

Then, d RTA and a RTA turn out to be both proportional 
to the non-dispersive intra-band conductivity cri ntra of 
Eq. which evaluates exactly to cri ntra = with 

Co = ^2 + 2fcTln(e“5fr + 1)]. 

Incidentally, we notice that the conductivity <j±, A , being 
markedly anisotropic, is proportional to ch 2 cri ntra . This 
makes its imaginary part change sign with increasing 
the frequency, in contrast to the purely inductive na¬ 
ture of <t±, A , whose imaginary part keeps a negative sign 
in the sampled frequency range. The plots of Fig. 0A) 
and 0B) illustrate such a behavior. Now, by definition, 
surface waves are sustained at the interface between ma¬ 
terials having different permittivities so that, in a trans¬ 
verse resonance circuit, their associated reactances cancel 
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FIG. 1. Real and imaginary parts of the scalar conductivities 
£Trta (Aj and Crta (B) calculated from the R!TA. model 0 - 
and reported vs the dimensionless frequency uj/qv ¥ . The two 
quantities are normalized to £o/gv F . 


out. Therefore, the contribution of the anisotropic term 
may strongly affect the polarization of surface plas- 
mons 22 of high rate of spatial variation. 

To include the diffusion currents, we introduce the fol¬ 
lowing tensor quantity 


s±(q) = 


zh. [ 

CO J (b — v±(k) • q 

1 st BZ 


(8) 


where 


f' 

J £ 


f 
J £ 


£±(k)-ii 


£±(k)-n 


I <Pkf' e 

1 st BZ 


(U A 4 


The BGK correction to the RTA conductivity can be put 
in the form 


°4gk(< 1) = [I + s±(q)] V^q). (9) 

This is equivalent to the result presented in Ref. and 
leads to correct the total conductivity as 

er(q) = <7 BG K(q) + ointerl- (10) 

The numerical results reported in the following are 
obtained by both a full-wave solver and a semi- 
analytical approach, as implemented respectively by the 
finite element method (FEM) and the method of mo¬ 
ments (MoM). The latter employs the usual free-space 
Green’s tensor G, in cylindrical coordinates r = (r, cp) 
and z, as the kernel of an integral operator relating the 
electric current density and field 20 21 \ 

E(r, z) = J d 2 r' j dz' G(r - r', * - *') • J(r', z'). 

Here, it should be noticed that the discretized currents in 
space are assumed as independent variables in the MoM 


calculation. Then, the solution for the resulting Elec¬ 
tric Field Integral Equation (EFIE) is generalized to the 
case of spatially dispersive material 16 * 20 * 21 *, by using the 
proper constitutive relation in the spatial domain, i.e., 

J s (r) = J d 2 r'u(r - r') • E s (r'). (11) 


In this equation, the planar current J s , and the tangent 
electric field E s , are sampled over the graphene surface. 

The non-local behavior of the current-field relation be¬ 
comes important when rapid field variations are involved. 
Transforming the conductivity in real space, the input 
spatial response of dispersive graphene is obtained as 


r ( r ) = I ^ 7^q°{q)M<i r ) 


■la; 


inter 


( 12 ) 
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where the scalar conductivities a and a include both the 
^-contributions from Eq. 0 or Eq. §. In Eq. @ 
the first and second terms contributing to the conduc¬ 
tance, are weighted by the zero-order and second-order 
Bessel functions, respectively. This means that the first 
addend of the conductivity is more sensitive to slow vary¬ 
ing fields with respect to the second addend. Equally im¬ 
portantly, the second term of the conductivity, containing 
off-diagonal matrix elements, is characterized by an an¬ 
gular dependence related to a current response, which is 
locally weighted by a “quadrupole” spatial distribution 
of the EM field. Consequently, unless very high angular 
and radial variations of the EM fields are concerned, the 
second term can be neglected. 

Let us now turn to the main application of the present 
work (Fig. [2]) that is a circular graphene disk of diameter 
D = 2/£, with a wire antenna of length L, placed just 
above its center. For this system, we can safely assume 
cylindrical symmetry. The gap between the antenna and 
the disk is L/200. In absence of angular variation of the 
EM excitation, the off diagonal terms vanish and no an¬ 
gular current arises. In addition, with a not too small 
tip-sample distance, and, thus, a not too strong EM field 


variation, the cr-term in Eq. (12) is expected to be dom¬ 
inant. Under these limiting conditions, Eq. © can be 
approximated by 


J r = J d 2 r'a(\r-r'\)E s (r') (13) 

where J r is the radial component of the surface current, 
and a(r) is the real space representation of a(g), which 
gives the (spatial) impulsive response of the current af¬ 
ter the EM excitation. Note that usual assumption of 
a thin hollow cylinder to approximate the wire antenna 
may affect the near field distribution between the tip 
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and graphene, but it does not limit the generality of the 
present analysis. 



FIG. 2. Wire antenna placed at sub-micrometric dis¬ 
tance (L/200) from a graphene patch 

For a fixed frequency, the effect of spatial dispersion 
increases with increasing the charge scattering-time. As 
a practical example, we take the cumulative integral of 
the impulsive response 

C(r) = J d 2 r'd(r') =2 tt J dr'r'd(r '), 

r'<r 

i.e., the current density response to a uniform unit elec¬ 
tric field within a circular area of radius r < R. Its 
profile, normalized to cri ntra , is shown in Fig. [3] for two 
different scattering times (r = 1, 2 ps), at / = 10 THz. 



FIG. 3. Cumulative integral of the impulsive response 
C/crintra as a function of radial position r, normalized to the 
free-space wavelength A; two different scattering times (r = 
1, 2 ps) are tested. 

In absence of dispersion, the response of Fig. [3] would 
be a real constant without spatial ripples. In presence 
of dispersion, the actual response is more complex, and 
some spatial ripples appear. These oscillations extend 
just to a fraction of the free-space wavelength A, which is 
typically comparable with the plasmon wavelength, and 
increases with increasing the charge lifetime. We expect 
that the above behavior of the system will reflect, numer¬ 
ically, on the solution of the EM field distribution. 


III. NUMERICAL RESULTS 


In the following, we present the solution for the EFIE 
directly in real space, focussing on both the full tensor 
form (12) and the scalar form (13) of the surface conduc¬ 
tivity. 


A. Non-Dispersive Analysis of Surface Plasmons 

We begin by considering an example of plasmon ex¬ 
citation without spatial dispersion. Let us define the 
graphene surface impedance Z s as the reciprocal of 
the non-dispersive Kubo conductivity a K introduced in 
Sec.|TTJ We then have Z s = l/a K or 1/Z S = a intra + cr int er, 
where the intra-band and inter-band terms have been re¬ 
spectively given in Eqs. © and ©• 


TIP 



FIG. 4. Plasmon distribution (normalized electric field) 
from a linear antenna coupled to a graphene disk (sketched 
in Fig. [2j vs the radial distance from the disk center and 
ReZ s /lmZ s . 

The plasmon distribution for a graphene disk of diam¬ 
eter D, coupled to an antenna of length L = D « c/2/, 
is reported in Fig. [4] vs the radial position within the disk 
and the inverse plasmon “ quality factor ”. The latter is 
defined as the ratio between imaginary and real parts of 
Z s calculated by the non-dispersive Kubo conductivity, 
i.e., lmZ s /ReZ s = — Rea K /Imcr K . In Fig. [4j both the op¬ 
eration frequency and the graphene chemical potential 
follow from the choice of the surface reactance ImZ s and 
plasmon quality-factor. For example, a reactance of 5 kU 
and a quality factor of 11 are associated to a chemical 
potential of 0.08 eV and a frequency of about 7.4 THz. 
We see that the spatial oscillations and propagation of 
the plasmon expire, at a progressively smaller radial dis¬ 
tance from the feeding tip, as ReZ s /ImZ s increases from 
0 to 0.15. 
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lm(Z s )=1.7 k Q 


lm(Z s )=2.0 k Q 


FIG. 5. (A) Input matching of the antenna vs surface re¬ 

actance; (B) iso-surface plot of the electric field (magnitude) 
radiated by the wire antenna, and plasmon excitation on the 
underlying graphene patch. 


Another interesting effect is pointed out in Fig. [5j 
where the overlap of the plasmon distribution with the 
near field under the tip determines the strength of the 
coupling and the amount of power transferred from the 
antenna to the plasmon. In order to emphasize this 
concept- and make the effect more evident- the losses are 
just neglected. More specifically, the surface resistivity 
is taken to be an independent variable, i.e., its value is 
not calculated from the Kubo conductivity as in Fig. [4] 
In addition, Z s is considered a purely imaginary quan¬ 
tity (no losses) and the reactance Im Z s is varied from 
about ~ 1.3 kf 1 to ~ 2.5 kf}. Thus, the geometric pa¬ 
rameters may be expressed in terms of the sizes of the 
antenna and the graphene disk relative to the vacuum 
wavelength A. In particular: (i) the diameter D of the 
patch and the length L of the wire antenna are related 
by D = L = 2.11c// = 2.11A; (ii) the air-gap distance 
d between the antenna and the graphene-patch is given 
by d = L/ 200; (iii) the radius a of the antenna is fixed 
to a = L/ 1000; (iv) a value of L/200 is chosen for the 
excitation-gap g in the middle of the antenna, where the 
voltage-source is applied; (v) The internal impedance Zq 


of the voltage applied to the excitation-gap of the an¬ 
tenna is set to 70 ft. 

It should be noted that the value of g can be scaled 
up or down provided that it remains much smaller than 
A and that the applied voltage is scaled reversely. In 
addition, the internal impedance value is equal to the 
real part of the input impedance Zll of the antenna 
in absence of graphene to have maximum power trans¬ 
fer. Fig. IA) shows that the antenna, which would 
work precisely at its resonance frequency in absence of 
graphene, can be more or less “detuned” by the graphene 
patch. The detuning depends on the graphene-antenna 
coupling: the higher the coupling strength, the higher 
the coupled reactive power, and the higher the reflection 
coefficient 511 at the input port of the antenna (indi¬ 
cated by arrows). 511 expresses the amount of power 
that is reflected back at the terminals of a voltage-source 
excitation, given by an infinitesimal electric dipole lo¬ 
cated at the center of the antenna. By definition, we 
have 511 = (Zll — Z 0 )/(Z11 + Z 0 ). The choice of the 
source impedance of the voltage excitation is arbitrary. 
However, as stated above, we have set Z 0 in order to have 
resonance in absence of graphene, and accordingly a bet¬ 
ter visualization of the perturbing effect of the graphene 
patch (see Sec. m Blfor more information and numerical 
values). Fig. §B) shows in details what is happening to 
the antenna near (right picture) and far from (left pic¬ 
ture) the detuning points, i.e., the arrows of Fig. §A). 
In particular, we have taken surface reactance values of 2 
and 1.7 kf}, corresponding respectively to input matching 
values of about —4 and — 13dB: in the former case, the 
radiation is drastically reduced by the plasmon coupling 
and a large amount of power is reflected back to the input 
port of the antenna. The results reported in Fig.[5|A) are 
derived from our MoM simulator 23 ^, whereas the plots of 
Fig. [5|B) are calculated by a full-wave EM solver (HFSS 
by Ansoft), which provides an independent validation of 
our implementation. 


B. Dispersive Characterization of Surface 
Plasmons 


In order to discuss the effect of dispersion on the 
strength of plasmon excitation, we focus on the two 
examples reported in Fig. [6j Differently from the re¬ 
sults shown in Fig. [5j here we take into account ohmic 
losses using the approximate complex conductivity of 


Eq. (13) within the RTA limit. Assuming a nominal 


frequency of 10 THz, we show the effect of spatial dis¬ 
persion on the input matching (Sll) of the antenna in 
resonant [L « A/2.11, Fig. |§A)] and non-resonant [L « 
A/12.5, Fig.[6](B)] conditions, respectively. All the reflec- 
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tion coefficients at the input terminal of the antenna are 
plotted as function of the chemical potential of graphene, 
i.e., the Fermi energy shift associated to the local doping 
level. 



FIG. 6. Reflection coefficient at the input terminal of the 
antenna with (solid lines) and without (dashed lines) spatial 
dispersion, for different relaxation times. The operating fre¬ 
quency is set to 10 THz. Two disk diameters are considered, 
namely D = L « A/2.11 (A), D = L « A/12.5 (B). 

Three typical values of charge-carrier relaxation times 
on the ps time-scale are tested (r = 0.5,1,2 ps). In¬ 
deed, the relaxation time in graphene may strongly de¬ 
pend on the quality of the sample related to the fabri¬ 
cation process. In particular, real samples are affected 
by the presence of grain boundaries, defects, multilayer 
regions, etc. For the above reason, some flexibility is 
needed in selecting the r parameter. As evident from 
Fig.§ the larger is the relaxation time the sharper are 
the Sll peaks that express the maxima and minima of the 
antenna-plasmon coupling. The antenna is fed by a volt¬ 
age source with an impedance equal to the real part of the 
input impedance of the antenna, which is about 70 O and 
2222 Q for the resonant and non-resonant cases, respec¬ 
tively. Clearly, the minima of reflection, corresponding 
to surface field resonances, have lower matching levels in 
the non-resonant case. In either resonant or non reso¬ 
nant condition, the dispersive reflection peaks are more 
broadened- and have lower maxima- with respect to the 


corresponding non-dispersive ones. These differences be¬ 
come more and more evident as the relaxation time in¬ 
creases, or, the quality of graphene gets better. 

In SMM, an electrostatic tuning of the charge density 
with a DC voltage applied to the microscope tipPcan pro¬ 
vide the local doping for near field applications. Fig. [7] 
shows the spatial distributions of the plasmonic wave on 
the graphene disk (for a normalized electric field) as func¬ 
tion of the radial position and the doping levels (values 
corresponding to those of Fig. |6|. The distributions are 
computed with and without spatial dispersion, assuming 
a relaxation time of 2 ps. Differences between the dis¬ 
persing and non dispersing curves are clearly observable, 
particularly in correspondence of the resonance peaks. 
Looking at the plots of Fig. [6] and Fig. [7| we see that 
the absolute effect of spatial dispersion is higher in the 
non-resonant case, where the size of the graphene disk is 
much smaller than the free-space wavelength A. 


non dispersive dispersive 



FIG. 7. Spatial distribution of the surface electric field 
E s (normalized to 1) as function of relative position (r/R) and 
chemical potential /i, with (B,D) and without (A,C) spatial 
dispersion effects, for f= 10 THz, and D = L ps A/2.11 (A,B), 
D = L?z A/12.5 (C,D) 


To go beyond the diagonal conductivity approximation 
of Eq. (13) in the RTA, we explicitly account for the an¬ 
gular variations of the EM field and, at the same time, 
we include the effect of diffusion currents. Thus, we next 
consider the full conductivity response of Eq. 
both the RTA and BGK approach. In Fig [8 


12 ] ) , within 
, we report 


the details of the input matching of the antenna in the 




































resonant case D « L = A/2.12, comparing the results 
obtained from the different levels of approximations dis¬ 
cussed in the present work. We find confirmation that the 
RTA conductivity of the graphene disk represents a sig¬ 
nificant improvement with respect to the non-dispersive 
conductivity of an infinite graphene sheet. In addition, 
we see that the contribution of the off diagonal tensor a 
leads to a non-negligible small correction to the diagonal 
tensor d, within the RTA. 



FIG. 8. Reflection coefficient (Sll) at input terminal of the 
antenna with the different approaches to dispersion discussed 
here. The following parameters are used: / = 10 THz, r = 
1 ps, and D — L « A/2.11. 


BGK form. 

IV. CONCLUSIONS 

We have used a semiclassical model derived from the 
Boltzmann transport equation to investigate the effect of 
spatial dispersion on the linear THz response of graphene, 
characterized by excitation of surface plasmons. The ex¬ 
citation source has been provided by the near field of an 
antenna radiating in proximity of a graphene micro-disk. 

We have characterized the role of spatial dispersion, 
with respect to a cylindrical system, within both the 
RTA and the BGK approach, obtaining meaningful and 
compact expressions of the full-tensor representing the 
constitutive relation of the graphene patch in real space. 

We have shown that the surface distribution of the 
field on the disk, and the macroscopic response of the 
antenna, is significantly affected by spatial dispersion in 
two distinct noteworthy examples, where the antenna has 
been set in resonant and non-resonant conditions. 

Although the role played by spatial dispersion was pre¬ 
viously clarified in propagation and modal analysis of 
plasmons 15 * 16 ^, here we have provided a focus on near field 
excitation problems, with potential fallout in important 
practical applications, concerning near field imaging. 


More importantly, we notice that the effect of electron 
diffusion included in the BGK model plays a significant 
role. Indeed the RTA and BGK expressions of a and a 
are remarkably different: dispersion effects appear to be 
under-estimated by the RTA conductivity [or even more 
simplified expressions as Eq. (13)] with respect to the 
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to Eq. [6j Now, consider the integral identity 


f 27T dip COS flip 
J 0 1 + Z COS -0 


27T(-l) n / VT^2 - 1 

VT^ V * 


which holds true for n = 0,1,... and Imz ^ 0, and can 
be derived from the databasd^. As special cases, we get 



dip 

1 — Z COS 0 


27r 


for 72 = 0 and 



dip cos 2ip 

1 — Z COS 0 



for 72 = 2. By a simple change of variable, we also have: 


APPENDIX: DISPERSIVE AND 
NON-DISPERSIVE RTA CONDUCTIVITY 


Using the polar coordinates k = k(cosd k , sind k ) and 
q = g(cos60,sin6^), the tensor product at the numerator 
of the BZ integral in Eq. 0 becomes 


v±(k)v±(k) = Vp 


cos 2 d k sin d k cos d k 
sin cos 7^ sin 2 d k 


l 

L 


dip cos 2 ip 


1 — z cos (ip — -0o) 
2?r dipsin2ip 


B(z) cos 2 * 00 , 
= B(z) sin2 , 0o. 


1 — zcos(ip — ipo) 

It follows that the dp -integral in Eq. (| 6 |) evaluates to 



dd k 


1 + 

cos 2d k sin 2d k 
sin 2 d k — cos 2 d k 

Co 

- v F gcos(i?fc - 0 q ) 


A(v F q/Cj) l 

Co 


u+ 

cos 2 d k 

sin 2d k 

\ . B(v F q/Q) 

cos 20 q sin 20 q 

r + 

sin 2d k 

— cos 2d k 

r + * 

sin 20 q — cos 20 q 


while the scalar products at the denominator reads: Plugging this in to Eq. we obtain Eq. ([t]) and the 

v±(k) • q = v F g cos{d k — 0 q ). Hence, Eq. [l] is turned scalar conductivities a^ a , a^ a . 




























